A meta-analysis of sex differences in animal personality: no evidence for greater male variability
Lauren Harrison, Daniel Noble & Michael Jennions
last updated: 16/09/2021
1 Supplementary Material - Data checks, analysis and meta-a models
2 Data Set-up
2.1 Load packages
# Clear work space
# rm(list=ls())
# Install CRAN packages
library("pacman")## Warning: package 'pacman' was built under R version 3.5.2# Install orchard plot and metaAidR packages from GitHub
devtools::install_github("itchyshin/orchard_plot", subdir = "orchaRd", force = TRUE, build_vignettes = TRUE); ## Error in get(genname, envir = envir) : object 'testthat_print' not found
##
##
checking for file ‘/private/var/folders/0b/pxghylq157gfhs1vrzdpx2gc0000gq/T/RtmpXG0rME/remotes6c851575cc49/itchyshin-orchard_plot-27d6281/orchaRd/DESCRIPTION’ ...
✓ checking for file ‘/private/var/folders/0b/pxghylq157gfhs1vrzdpx2gc0000gq/T/RtmpXG0rME/remotes6c851575cc49/itchyshin-orchard_plot-27d6281/orchaRd/DESCRIPTION’
##
─ preparing ‘orchaRd’:
##
checking DESCRIPTION meta-information ...
✓ checking DESCRIPTION meta-information
##
─ installing the package to build vignettes
##
creating vignettes ...
✓ creating vignettes (1.4s)
##
─ checking for LF line-endings in source and make files and shell scripts
##
─ checking for empty or unneeded directories
## ─ looking to see if a ‘data/datalist’ file should be added
##
─ building ‘orchaRd_0.0.0.9000.tar.gz’
##
## devtools::install_github("daniel1noble/metaAidR");
pacman::p_load(knitr, metafor, dplyr, kableExtra, tidyverse, rotl, phytools, GGally, R.rsp, patchwork, devtools, robumeta, ape, geiger, phytools, phangorn, rlist, orchaRd, metaAidR, corrplot, stringr)
# set working directory
# setwd("~/Documents/GitHub/sex_meta/")
# Source our own functions
source("./R/func.R")
# Set the rerun object to FALSE so that you don't need to re-run all models again. Some take quite a lot of time to run. If FALSE, it will just re-load saved output.
rerun_models = FALSE3 Organising data for analysis
## load our original pers dataset and our dataset with all of our sexual selection info
pers <- read.csv("./data/pers_data.csv", stringsAsFactors = FALSE)
bodysize <- read.csv("./data/bodysize_SSD.csv", stringsAsFactors = FALSE)
## Merge the two by spp_names columns
pers <- merge(x = pers,
y = bodysize[,c("species_name", "SSD_index", "mating_system")],
by="species_name", all.x=TRUE, no.dups = TRUE)
## Select only the relevant columns to make life easier
pers_new <- pers %>%
select(study_ID, year, species_name, SSD_index, taxo_group, data_type, personality_trait, male_n, male_mean_conv,
male_sd_conv, female_n, female_mean_conv, female_sd_conv, depend, directionality, spp_name_phylo, mating_system,
age, population, study_environment, study_type, measurement_type)
## Add in observation level random effect (metafor doesn't do this, need to do it manually)
pers_new <- pers_new %>%
group_by(taxo_group) %>%
mutate(obs = 1:length(study_ID)) 4 Calculating effect sizes
4.1 means - SMD using Hedges' g and variability - lnCVR
## SMD (Hedges' g)
pers_new <- escalc(measure = "SMD",
n1i = male_n, n2i = female_n,
m1i = male_mean_conv, m2i = female_mean_conv,
sd1i = male_sd_conv, sd2i = female_sd_conv, data = pers_new, var.names=c("SMD_yi","SMD_vi"), append = TRUE)
## lnCVR
pers_new <- escalc(measure = "CVR",
n2i = female_n, n1i = male_n,
m2i = female_mean_conv, m1i = male_mean_conv,
sd2i = female_sd_conv, sd1i = male_sd_conv, data = pers_new, var.names=c("CVR_yi","CVR_vi"))
# we have some NAs where one or both sexes have a value of 0 for either mean or SD. Will be easiest to just remove these.
# Exclude NAs
pers_new <- pers_new %>%
filter(!is.na(CVR_yi), !is.na(SMD_yi))
dim(pers_new) # check they've been removed with no issues5 Data Checks
5.1 mean-variance relationship in our dataset
Looking at the strength of the correlation between the mean and SD to check that using lnCVR as a measure of variability is valid.
If the mean and SD are NOT strongly correlated then using lnCVR is pointless.
# females and males seperately because they are in different columns
# use ggplot to make a scatterplot of females
fem <- ggplot(pers_new, aes(x = female_mean_conv, y = female_sd_conv)) + geom_point()
# on log scale
fem + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10')
# mean and SD on log scale to calculate correlation
logfemale_mean <- log(pers_new$female_mean_conv)
logfemale_SD <- log(pers_new$female_sd_conv)
# correlation between mean and SD
cor(logfemale_mean, logfemale_SD) #0.91## [1] 0.9182668 # Males
# use ggplot to make a scatterplot of females
male <- ggplot(pers_new, aes(x = male_mean_conv, y = male_sd_conv)) + geom_point()
# on log scale
male + scale_x_continuous(trans = 'log10') + scale_y_continuous(trans = 'log10')
# mean and SD on log scale to calculate correlation
logmale_mean <- log(pers_new$male_mean_conv)
logmale_SD <- log(pers_new$male_sd_conv)
# correlation between mean and SD
cor(logmale_mean, logmale_SD) #0.90## [1] 0.90712255.2 Checking for outliers and removing weird effect sizes
This is an important data checking step - here we can identify whether data has been entered or reported incorrectly (i.e. outliers)
First, let's look at the funnel plots for lnCVR and SMD. NOTE: these funnels have had our 2 big outliers already removed.
#funnel plot for lnCVR
funnel(x = pers_new$CVR_yi, vi = pers_new$CVR_vi, yaxis="seinv", xlim = c(-10, 10))
SMD
#funnel plot for SMD
funnel(x = pers_new$SMD_yi, vi = pers_new$SMD_vi, yaxis="seinv")
text(as.character(pers_new$study_ID), x = pers_new$SMD_yi, y = 1/sqrt(pers_new$SMD_vi))
Removing outliers
pers_new %>%
filter(study_ID == "P015" & SMD_yi < -15) # P015 has 1 large effect size, remove?
# filter out that large effect size
pers_new <- pers_new %>%
filter(!study_ID == "P015" | !obs == "509")
dim(pers_new)
# checking SMD outliers - inverse SE > 14
funnel(x = pers_new$SMD_yi, vi = pers_new$SMD_vi, yaxis="seinv")
text(as.character(pers_new$obs), x = pers_new$SMD_yi, y = 1/sqrt(pers_new$SMD_vi), offset = 0.8)
# checking lnCVR outliers
funnel(x = pers_new$CVR_yi, vi = pers_new$CVR_vi, yaxis="seinv", xlim = c(-10, 10))
text(as.character(pers_new$obs), x = pers_new$CVR_yi, y = 1/sqrt(pers_new$CVR_vi), offset = 0.8)
# Some measures are more physiological/not personality than personality, so probably wise to remove these before we run the models:
# P029 - obs 22, 23, 32
# P084 - obs 59, 62, 63, 65, 68, 70, 71, 72, 74
# P060 - obs 216, 217
# P211 - obs 230, 245
# P117 - obs 397, 393, 402, 414
# P197 - obs 541, 544, 546, 547
# P069 - obs 669, 672, 673, 682, 683, 684, 686, 694
# remove these by study
pers_new <- pers_new %>% filter(!study_ID == "P029" | !obs %in% c("21", "25", "26", "28", "31"))
pers_new <- pers_new %>% filter(!study_ID == "P084" | !obs %in% c("59", "62", "63", "65", "68", "70", "71", "72", "74"))
pers_new <- pers_new %>% filter(!study_ID == "P060" | !obs %in% c("216", "217"))
pers_new <- pers_new %>% filter(!study_ID == "P211" | !obs %in% c("230", "245"))
pers_new <- pers_new %>% filter(!study_ID == "P117" | !obs %in% c("397", "393", "402", "414"))
pers_new <- pers_new %>% filter(!study_ID == "P197" | !obs %in% c("541", "544", "546", "547"))
pers_new <- pers_new %>% filter(!study_ID == "P197" | !obs %in% c("669", "672", "673", "682", "683", "684", "686", "694"))
pers_new <- pers_new %>% filter(!study_ID == "P041" | !obs %in% c("120", "124"))
dim(pers_new) # check they've been removed without issue
# after checking the data, there are a few effect sizes that might be driving weird results so let's drop them and see
pers_new <- pers_new %>% filter(!study_ID == "P100" | !obs == "519") # big outlier5.3 Flip signs of effects for SMD
The directional meaning of effect sizes vary depending on the specific units and trait being measured. The data has a directionality column that tells one if the meaning should be reversed (1) or left the same.
pers_new$directionality <- ifelse(is.na(pers_new$directionality), 0, 1)
pers_new$SMD_yi_flip <- ifelse(pers_new$directionality == 1, pers_new$SMD_yi*(-1), pers_new$SMD_yi)6 Prepare the phylogenetic trees
We constructed seperate phylogenetic trees for each taxonomic group. The tree for birds was constructed using BirdTree.org, the rest were constructed using TimeTree.org. We'll use these trees for multi-level meta-analytic models throughout the analysis.
# Find all tree file names
tree_files <- paste0("./trees/", list.files("./trees"))[-1]
# Bird tree has been constructed already, just need to get trees for the rest of the taxo groups
trees <- lapply(tree_files, function(x) read.tree(x))
names <- gsub("~./trees/", "", tree_files)
names(trees) <- names
# Plot the trees and see how they look
par(mfrow = c(1,5), mar = c(1,1,1,1))
lapply(trees, function(x) plot(x, cex = 1))
## $`./trees/bird_species.nwk`
## $`./trees/bird_species.nwk`$type
## [1] "phylogram"
##
## $`./trees/bird_species.nwk`$use.edge.length
## [1] TRUE
##
## $`./trees/bird_species.nwk`$node.pos
## [1] 1
##
## $`./trees/bird_species.nwk`$node.depth
## [1] 1
##
## $`./trees/bird_species.nwk`$show.tip.label
## [1] TRUE
##
## $`./trees/bird_species.nwk`$show.node.label
## [1] FALSE
##
## $`./trees/bird_species.nwk`$font
## [1] 3
##
## $`./trees/bird_species.nwk`$cex
## [1] 1
##
## $`./trees/bird_species.nwk`$adj
## [1] 0
##
## $`./trees/bird_species.nwk`$srt
## [1] 0
##
## $`./trees/bird_species.nwk`$no.margin
## [1] FALSE
##
## $`./trees/bird_species.nwk`$label.offset
## [1] 0
##
## $`./trees/bird_species.nwk`$x.lim
## [1] 0.000 1716.674
##
## $`./trees/bird_species.nwk`$y.lim
## [1] 1 109
##
## $`./trees/bird_species.nwk`$direction
## [1] "rightwards"
##
## $`./trees/bird_species.nwk`$tip.color
## [1] "black"
##
## $`./trees/bird_species.nwk`$Ntip
## [1] 109
##
## $`./trees/bird_species.nwk`$Nnode
## [1] 108
##
## $`./trees/bird_species.nwk`$root.time
## NULL
##
## $`./trees/bird_species.nwk`$align.tip.label
## [1] FALSE
##
##
## $`./trees/fish_species.nwk`
## $`./trees/fish_species.nwk`$type
## [1] "phylogram"
##
## $`./trees/fish_species.nwk`$use.edge.length
## [1] TRUE
##
## $`./trees/fish_species.nwk`$node.pos
## [1] 1
##
## $`./trees/fish_species.nwk`$node.depth
## [1] 1
##
## $`./trees/fish_species.nwk`$show.tip.label
## [1] TRUE
##
## $`./trees/fish_species.nwk`$show.node.label
## [1] FALSE
##
## $`./trees/fish_species.nwk`$font
## [1] 3
##
## $`./trees/fish_species.nwk`$cex
## [1] 1
##
## $`./trees/fish_species.nwk`$adj
## [1] 0
##
## $`./trees/fish_species.nwk`$srt
## [1] 0
##
## $`./trees/fish_species.nwk`$no.margin
## [1] FALSE
##
## $`./trees/fish_species.nwk`$label.offset
## [1] 0
##
## $`./trees/fish_species.nwk`$x.lim
## [1] 0.00 35169.86
##
## $`./trees/fish_species.nwk`$y.lim
## [1] 1 22
##
## $`./trees/fish_species.nwk`$direction
## [1] "rightwards"
##
## $`./trees/fish_species.nwk`$tip.color
## [1] "black"
##
## $`./trees/fish_species.nwk`$Ntip
## [1] 22
##
## $`./trees/fish_species.nwk`$Nnode
## [1] 21
##
## $`./trees/fish_species.nwk`$root.time
## NULL
##
## $`./trees/fish_species.nwk`$align.tip.label
## [1] FALSE
##
##
## $`./trees/invert_species.nwk`
## $`./trees/invert_species.nwk`$type
## [1] "phylogram"
##
## $`./trees/invert_species.nwk`$use.edge.length
## [1] TRUE
##
## $`./trees/invert_species.nwk`$node.pos
## [1] 1
##
## $`./trees/invert_species.nwk`$node.depth
## [1] 1
##
## $`./trees/invert_species.nwk`$show.tip.label
## [1] TRUE
##
## $`./trees/invert_species.nwk`$show.node.label
## [1] FALSE
##
## $`./trees/invert_species.nwk`$font
## [1] 3
##
## $`./trees/invert_species.nwk`$cex
## [1] 1
##
## $`./trees/invert_species.nwk`$adj
## [1] 0
##
## $`./trees/invert_species.nwk`$srt
## [1] 0
##
## $`./trees/invert_species.nwk`$no.margin
## [1] FALSE
##
## $`./trees/invert_species.nwk`$label.offset
## [1] 0
##
## $`./trees/invert_species.nwk`$x.lim
## [1] 0.000 8492.817
##
## $`./trees/invert_species.nwk`$y.lim
## [1] 1 42
##
## $`./trees/invert_species.nwk`$direction
## [1] "rightwards"
##
## $`./trees/invert_species.nwk`$tip.color
## [1] "black"
##
## $`./trees/invert_species.nwk`$Ntip
## [1] 42
##
## $`./trees/invert_species.nwk`$Nnode
## [1] 41
##
## $`./trees/invert_species.nwk`$root.time
## NULL
##
## $`./trees/invert_species.nwk`$align.tip.label
## [1] FALSE
##
##
## $`./trees/mammal_species.nwk`
## $`./trees/mammal_species.nwk`$type
## [1] "phylogram"
##
## $`./trees/mammal_species.nwk`$use.edge.length
## [1] TRUE
##
## $`./trees/mammal_species.nwk`$node.pos
## [1] 1
##
## $`./trees/mammal_species.nwk`$node.depth
## [1] 1
##
## $`./trees/mammal_species.nwk`$show.tip.label
## [1] TRUE
##
## $`./trees/mammal_species.nwk`$show.node.label
## [1] FALSE
##
## $`./trees/mammal_species.nwk`$font
## [1] 3
##
## $`./trees/mammal_species.nwk`$cex
## [1] 1
##
## $`./trees/mammal_species.nwk`$adj
## [1] 0
##
## $`./trees/mammal_species.nwk`$srt
## [1] 0
##
## $`./trees/mammal_species.nwk`$no.margin
## [1] FALSE
##
## $`./trees/mammal_species.nwk`$label.offset
## [1] 0
##
## $`./trees/mammal_species.nwk`$x.lim
## [1] 0.000 1711.683
##
## $`./trees/mammal_species.nwk`$y.lim
## [1] 1 46
##
## $`./trees/mammal_species.nwk`$direction
## [1] "rightwards"
##
## $`./trees/mammal_species.nwk`$tip.color
## [1] "black"
##
## $`./trees/mammal_species.nwk`$Ntip
## [1] 46
##
## $`./trees/mammal_species.nwk`$Nnode
## [1] 45
##
## $`./trees/mammal_species.nwk`$root.time
## NULL
##
## $`./trees/mammal_species.nwk`$align.tip.label
## [1] FALSE
##
##
## $`./trees/reptile_species.nwk`
## $`./trees/reptile_species.nwk`$type
## [1] "phylogram"
##
## $`./trees/reptile_species.nwk`$use.edge.length
## [1] TRUE
##
## $`./trees/reptile_species.nwk`$node.pos
## [1] 1
##
## $`./trees/reptile_species.nwk`$node.depth
## [1] 1
##
## $`./trees/reptile_species.nwk`$show.tip.label
## [1] TRUE
##
## $`./trees/reptile_species.nwk`$show.node.label
## [1] FALSE
##
## $`./trees/reptile_species.nwk`$font
## [1] 3
##
## $`./trees/reptile_species.nwk`$cex
## [1] 1
##
## $`./trees/reptile_species.nwk`$adj
## [1] 0
##
## $`./trees/reptile_species.nwk`$srt
## [1] 0
##
## $`./trees/reptile_species.nwk`$no.margin
## [1] FALSE
##
## $`./trees/reptile_species.nwk`$label.offset
## [1] 0
##
## $`./trees/reptile_species.nwk`$x.lim
## [1] 0.000 2656.912
##
## $`./trees/reptile_species.nwk`$y.lim
## [1] 1 10
##
## $`./trees/reptile_species.nwk`$direction
## [1] "rightwards"
##
## $`./trees/reptile_species.nwk`$tip.color
## [1] "black"
##
## $`./trees/reptile_species.nwk`$Ntip
## [1] 10
##
## $`./trees/reptile_species.nwk`$Nnode
## [1] 9
##
## $`./trees/reptile_species.nwk`$root.time
## NULL
##
## $`./trees/reptile_species.nwk`$align.tip.label
## [1] FALSE# Checking trees to ensure we only include species in the current dataset
# Check that they are ultrametric
lapply(trees, function(x) is.ultrametric(x))
# Check that all names in the phylogeny are also in the data
taxa_data_list <- split(pers_new, pers_new$taxo_group)
other_groups <- mapply(x = taxa_data_list,
y = trees,
function(x,y) tree_checks(x,y, "spp_name_phylo", type = "checks"))
# Print out each taxon group
for(i in colnames(other_groups)){
print(i)
print(other_groups[,i] )
}
# Now to prune trees so that we get tree names that match with species in data
pruned_trees <- mapply(x = taxa_data_list,
y = trees,
function(x,y) tree_checks(x,y, "spp_name_phylo", type = "prune"))
# Check that this has been done correctly
re_checks <- mapply(x = taxa_data_list,
y = pruned_trees,
function(x,y) tree_checks(x,y, "spp_name_phylo", type = "checks"))
for(i in colnames(re_checks)){
print(i)
print(re_checks[,i] )
}
# Extract the phylogenetic correlation matrices
phylo_vcv <- lapply(pruned_trees, function(x) vcv(x, corr = TRUE)) # these matrices are used in the meta-a models7 Sensitivity analysis - checking score data
Before we begin, we need to run a sensitivity analysis to see if score data is ok to use. With these models, we are just including score as a moderator term to compare with the rest of the dataset (some of which has already been transformed, we just can't do that with scores). Model summaries are also presented in Supplementary Table S2.
Our score sensitivity model:
# model:
sensitivity_mod1_score <- meta_model_fits(pers_new, phylo_vcv, type = "score")
# Extract the SMD and lnCVR results
smd_mods_score <- sensitivity_mod1_score["SMD",]
lnCVR_mods_score <- sensitivity_mod1_score["lnCVR",] # inverts significant
# Because invert score data is significantly different, we need to remove these effect sizes before running our models
# filter out invert scores from pers dataset
pers_new <- pers_new %>%
filter(score != "score" | taxo_group != "invertebrate")
dim(pers_new)
# How many effect sizes, unique studies and different species are we left with?
data.frame(pers_new %>%
summarise(n = n(), studies = length(unique(study_ID)), species = length(unique(spp_name_phylo))))8 Meta-Analysis Models
Let's run the first bunch of models on the whole dataset. We'll start off with intercept-only multi-level meta-analytic models, then move to multi-level meta-regression models (personality traits, and SSD). The functions in func.R should be consulted to see precisely what models are being fit across the taxonomic groups.
8.1 Intercept-only MLMA models
Complete model summaries are also presented in Supplementary Table S14.
# First we will fit our MLMA intercept only models, across each taxo group.
# we can use this function to just read the saved model output instead of re-running the model, which takes a while
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMA_models <- meta_model_fits(pers_new, phylo_vcv, type = "int")
saveRDS(MLMA_models, "./output/MLMA_models_int")
}else{
MLMA_models <- readRDS("./output/MLMA_models_int")
}
# View model results
split_taxa <- split(pers_new, pers_new$taxo_group)
smd_mods <- MLMA_models["SMD",]
lnCVR_mods <- MLMA_models["lnCVR",]8.2 I2 estimates of heterogeneity - intercept models
Study_ID is the between study heterogeneity, Phylo tells us if there is a phylogenetic signal and the strength of that signal. Total I2 is testing how much heterogeneity we have beyond sampling variance
# From these models we can get I2 estimates:
birds_smd = I2(smd_mods[[1]], v = split_taxa[[1]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
birds_CVR = I2(lnCVR_mods[[1]], v = split_taxa[[1]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
fish_smd = I2(smd_mods[[2]], v = split_taxa[[2]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
fish_CVR = I2(lnCVR_mods[[2]], v = split_taxa[[2]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
invert_smd = I2(smd_mods[[3]], v = split_taxa[[3]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
invert_CVR = I2(lnCVR_mods[[3]], v = split_taxa[[3]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
mammal_smd = I2(smd_mods[[4]], v = split_taxa[[4]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
mammal_CVR = I2(lnCVR_mods[[4]], v = split_taxa[[4]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
reptile_smd = I2(smd_mods[[5]], v = split_taxa[[5]]$SMD_vi, phylo = "spp_name_phylo", obs = "obs")
reptile_CVR = I2(lnCVR_mods[[5]], v = split_taxa[[5]]$CVR_vi, phylo = "spp_name_phylo", obs = "obs")
# Now that we have our list of models, we can extract the estimates, CIs and prediction intervals
MLMA_estimates_SMD <- plyr::ldply(lapply(smd_mods, function(x) print(mod_results(x, mod = "Int"))))
MLMA_estimates_SMD MLMA_estimates_lnCVR <- plyr::ldply(lapply(lnCVR_mods, function(x) print(mod_results(x, mod = "Int"))))
MLMA_estimates_lnCVRExtract p-values from these models to use later when adjusting them for multiple testing:
# taking p-values from models for False Discovery Rate p-value adjustment
p.SMD_intercept <- unlist(lapply(smd_mods, function(x) x$pval))
p.SMD_intercept## bird fish invertebrate mammal reptilia
## 0.3560563 0.3900624 0.1778479 0.7143065 0.4793203 p.lnCVR_intercept <- unlist(lapply(lnCVR_mods, function(x) x$pval))
p.lnCVR_intercept## bird fish invertebrate mammal reptilia
## 0.5812731 0.9117446 0.8659517 0.6776659 0.35541548.3 Personality trait MLMR models
These models include personality trait type as a moderator. Please note that we estimate the mean for each of the categorical levels because we are not really interested in whether the means differ, but whether or not males and females differ in any of these traits.
Complete model summaries are presented in Supplementary Table S15.
# we can just reload saved model outputs here to save time
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMR_models_pers_trait <- meta_model_fits(pers_new, phylo_vcv, type = "pers")
saveRDS(MLMR_models_pers_trait, "./output/MLMR_models_pers_trait")
} else{
MLMR_models_pers_trait <- readRDS("./output/MLMR_models_pers_trait")
}
# Extract the SMD and lnCVR results
smd_mods_pers <- MLMR_models_pers_trait["SMD",]
lnCVR_mods_pers <- MLMR_models_pers_trait["lnCVR",]
# these model objects are used to make the orchard plots shown in Figures 2-6Get prediction intervals for personality trait models:
# Get the combined estimates from them all
MLMA_estimates_SMD_pers <- plyr::ldply(lapply(smd_mods_pers, function(x) print(mod_results(x, mod = "personality_trait"))))
MLMA_estimates_SMD_pers MLMA_estimates_lnCVR_pers <- plyr::ldply(lapply(lnCVR_mods_pers, function(x) print(mod_results(x, mod = "personality_trait"))))
MLMA_estimates_lnCVR_pers # Add in n and k to these dataframes
n_k<- pers_new %>%
group_by(taxo_group, personality_trait) %>%
summarise(n = n(), spp = length(unique(spp_name_phylo)), k = length(unique(study_ID)))
# Summary of model estimates with number of studies, species and effect sizes included
MLMA_estimates_SMD_pers <- data.frame(MLMA_estimates_SMD_pers, n_k[,c("n", "spp", "k")])
MLMA_estimates_SMD_pers MLMA_estimates_lnCVR_pers <- data.frame(MLMA_estimates_lnCVR_pers, n_k[,c("n", "spp", "k")])
MLMA_estimates_lnCVR_persExtract p-values from models for multiple testing adjustment later:
# extract p-values for multiple testing
p.SMD_pers <- unlist(lapply(smd_mods_pers, function(x) x$pval))
p.SMD_pers## bird1 bird2 bird3 bird4 bird5
## 0.57653039 0.30479363 0.27751708 0.65883812 0.01348632
## fish1 fish2 fish3 fish4 fish5
## 0.51728620 0.27173344 0.70897900 0.26038453 0.54643814
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5
## 0.21528997 0.23621794 0.22305017 0.86210097 0.20535070
## mammal1 mammal2 mammal3 mammal4 mammal5
## 0.42949385 0.64539838 0.47941582 0.89056692 0.74385107
## reptilia1 reptilia2 reptilia3 reptilia4 reptilia5
## 0.83127102 0.45179068 0.38700906 0.03843513 0.98694696 p.lnCVR_pers <- unlist(lapply(lnCVR_mods_pers, function(x) x$pval))
p.lnCVR_pers## bird1 bird2 bird3 bird4 bird5
## 0.66827246 0.63819433 0.71794905 0.02062033 0.56829592
## fish1 fish2 fish3 fish4 fish5
## 0.89295129 0.08859994 0.67087141 0.63106746 0.02038126
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5
## 0.54416559 0.80512014 0.72746311 0.21252880 0.56927261
## mammal1 mammal2 mammal3 mammal4 mammal5
## 0.46265634 0.45330729 0.82211806 0.72515472 0.76995584
## reptilia1 reptilia2 reptilia3 reptilia4 reptilia5
## 0.49192601 0.06178802 0.20080617 0.25517952 0.909719838.4 Personality trait x SSD MLMR models
Now let's look at how SSD interacts with personality trait type. Here we are not estimating an intercept either, so each intercept varies by trait category and each slope as well. Note that there are lots of warnings, but these are the result of many levels not being present in taxa groups. We chose not to scale SSD_index because it is easier (and biologically relevant) to interpret SSD when it is 0 (when males and females are the same size), and when SSD is positive (when males are larger than females).
Model summaries are presented in Supplementary Table S17.
# again, we can just reload our saved model output here
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMR_models_pers_SSD <- meta_model_fits(pers_new, phylo_vcv, type = "pers_SSD")
saveRDS(MLMR_models_pers_SSD, "./output/MLMR_models_pers_SSD")
} else{
MLMR_models_pers_SSD <- readRDS("./output/MLMR_models_pers_SSD")
}
# Extract the SMD and lnCVR results
smd_mods_pers_SSD <- MLMR_models_pers_SSD["SMD",]
lnCVR_mods_pers_SSD <- MLMR_models_pers_SSD["lnCVR",]Get the prediction intervals for our interaction models:
# extract estimates using modified function in func.R file:
# SMD
MLMA_estimates_SMD_SSD <- plyr::ldply(lapply(smd_mods_pers_SSD, function(x)
print(mod_results_new(x, mod_cat = "personality_trait", mod_cont = "SSD_index", type = "zero"))))
MLMA_estimates_SMD_SSD # lnCVR
MLMA_estimates_lnCVR_SSD <- plyr::ldply(lapply(lnCVR_mods_pers_SSD, function(x)
print(mod_results_new(x, mod_cat = "personality_trait", mod_cont = "SSD_index", type = "zero"))))
MLMA_estimates_lnCVR_SSD # Table to get species numbers, no. studies and no. effect sizes:
data.frame(pers_new %>%
group_by(taxo_group, personality_trait) %>%
filter(!is.na(SSD_index))%>%
summarise(n = n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))8.5 SSD subset models
Because we aren't really interested in how each trait type differs from each other, we need to run our SSD models on subsets of the data where we can get the mean estimates for individual trait types and for SSD. Model summaries are presented in Supplementary Table S16.
NOTE: Since we are conducting our meta-regression at the species level (the level at which we can assume effect sizes are independent), any personality trait with fewer than 10 species needs to be dropped to look at interactions between SSD and personality. Having a minimum of 10 studies etc. is the rule of thumb for meta-regressions (e.g. see Borenstein et al Intro to Meta-A)
8.5.1 NEED TO DROP:
- ALL REPTILES - GONE, NOT ENOUGH SPECIES
- BIRDS - EVERYTHING BUT BOLDNESS
- FISH - ACTIVITY, EXPLORATION & SOCIALITY
- INVERTS - SOCIALITY, EXPLORATION & AGGRESSION
- MAMMALS - SOCIALITY
8.5.2 Mammals
SSD for activity, boldness, aggression and exploration.
# 1. MAMMALS
# First, we need to subset our pers dataset by taxo group to drop the unwanted levels.
# a. activity
pers_new_mammal_activity <- as.data.frame(pers_new %>%
filter(personality_trait == "activity") %>%
filter(taxo_group == "mammal"))
# b. boldness
pers_new_mammal_boldness <- as.data.frame(pers_new %>%
filter(personality_trait == "boldness") %>%
filter(taxo_group == "mammal"))
# c. aggression
pers_new_mammal_aggression <- as.data.frame(pers_new %>%
filter(personality_trait == "aggression") %>%
filter(taxo_group == "mammal"))
# d. exploration
pers_new_mammal_exploration <- as.data.frame(pers_new %>%
filter(personality_trait == "exploration") %>%
filter(taxo_group == "mammal"))
# Extract the phylogenetic correlation matrices
phylo_vcv_mammal <- phylo_vcv[[4]] Activity:
# a. activity
#SMD
MLMR_mods_pers_SSD_mammal_activity_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_activity)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data. MLMR_mods_pers_SSD_mammal_activity_SMD##
## Multivariate Meta-Analysis Model (k = 83; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1000 0.3163 14 no study_ID no
## sigma^2.2 2.8075 1.6755 12 no spp_name_phylo yes
## sigma^2.3 0.1914 0.4375 83 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 320.2151, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 81) = 5.0697, p-val = 0.0271
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.5040 1.2531 0.4022 81 0.6886 -1.9892 2.9972
## SSD_index -2.2054 0.9795 -2.2516 81 0.0271 -4.1543 -0.2565 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_SSD_mammal_activity_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_activity)## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data. MLMR_mods_pers_SSD_mammal_activity_lncvr##
## Multivariate Meta-Analysis Model (k = 83; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0293 0.1710 14 no study_ID no
## sigma^2.2 0.0001 0.0079 12 no spp_name_phylo yes
## sigma^2.3 0.0607 0.2465 83 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 81) = 146.2596, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 81) = 0.1324, p-val = 0.7169
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0518 0.0988 0.5236 81 0.6020 -0.1449 0.2484
## SSD_index 0.1248 0.3430 0.3638 81 0.7169 -0.5577 0.8073
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Boldness:
# b. boldness
#SMD
MLMR_mods_pers_SSD_mammal_bold_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_boldness)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.## Warning: Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_mammal_bold_SMD ##
## Multivariate Meta-Analysis Model (k = 163; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0088 0.0938 26 no study_ID no
## sigma^2.2 0.0000 0.0050 26 no spp_name_phylo yes
## sigma^2.3 0.1707 0.4132 163 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 161) = 405.7659, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 161) = 1.3101, p-val = 0.2541
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0739 0.0774 0.9547 161 0.3412 -0.0789 0.2267
## SSD_index -0.1686 0.1473 -1.1446 161 0.2541 -0.4594 0.1223
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_SSD_mammal_bold_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_boldness)## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_mammal_bold_lncvr##
## Multivariate Meta-Analysis Model (k = 163; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0029 0.0543 26 no study_ID no
## sigma^2.2 0.0000 0.0027 26 no spp_name_phylo yes
## sigma^2.3 0.0211 0.1452 163 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 161) = 177.7499, p-val = 0.1737
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 161) = 2.1066, p-val = 0.1486
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0114 0.0513 0.2230 161 0.8238 -0.0899 0.1128
## SSD_index 0.1316 0.0907 1.4514 161 0.1486 -0.0475 0.3106
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Aggression:
# c. aggression
#SMD
MLMR_mods_pers_SSD_mammal_aggression_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_aggression)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.## Warning: Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_mammal_aggression_SMD##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0000 0.0039 15 no study_ID no
## sigma^2.2 0.6852 0.8277 13 no spp_name_phylo yes
## sigma^2.3 0.1430 0.3781 85 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 312.3189, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 83) = 4.2481, p-val = 0.0424
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.1036 0.6014 -0.1722 83 0.8637 -1.2997 1.0926
## SSD_index 1.4134 0.6858 2.0611 83 0.0424 0.0495 2.7774 *
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_SSD_mammal_aggression_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_aggression)## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_mammal_aggression_lncvr##
## Multivariate Meta-Analysis Model (k = 85; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1790 0.4231 15 no study_ID no
## sigma^2.2 0.0000 0.0052 13 no spp_name_phylo yes
## sigma^2.3 0.1539 0.3922 85 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 83) = 202.3514, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 83) = 0.0118, p-val = 0.9138
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0992 0.1534 0.6467 83 0.5196 -0.2059 0.4042
## SSD_index -0.0756 0.6959 -0.1086 83 0.9138 -1.4596 1.3084
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Exploration:
# d. exploration
#SMD
MLMR_mods_pers_SSD_mammal_explore_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_exploration)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data. MLMR_mods_pers_SSD_mammal_explore_SMD##
## Multivariate Meta-Analysis Model (k = 213; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0504 0.2244 19 no study_ID no
## sigma^2.2 0.0000 0.0048 16 no spp_name_phylo yes
## sigma^2.3 0.1331 0.3649 213 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 211) = 658.4587, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 211) = 0.0351, p-val = 0.8516
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0016 0.0914 -0.0173 211 0.9862 -0.1817 0.1786
## SSD_index -0.0522 0.2786 -0.1873 211 0.8516 -0.6015 0.4971
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_SSD_mammal_explore_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal_exploration)## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data. MLMR_mods_pers_SSD_mammal_explore_lncvr##
## Multivariate Meta-Analysis Model (k = 213; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0198 0.1406 19 no study_ID no
## sigma^2.2 0.0265 0.1628 16 no spp_name_phylo yes
## sigma^2.3 0.0323 0.1799 213 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 211) = 361.1620, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 211) = 0.2658, p-val = 0.6067
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0595 0.1507 -0.3951 211 0.6932 -0.3566 0.2375
## SSD_index 0.1324 0.2567 0.5156 211 0.6067 -0.3737 0.6384
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.5.3 Birds
SSD for boldness only.
# 2. BIRDS
# subset dataset
pers_new_bird <- as.data.frame(pers_new %>%
filter(personality_trait == "boldness" & taxo_group == "bird"))
# phylo_vcv birds only
phylo_vcv_bird <- phylo_vcv[[1]]Boldness:
# SMD
MLMR_mods_pers_SSD_bird_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"),
test = "t", data = pers_new_bird)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.## Warning: Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_bird_SMD##
## Multivariate Meta-Analysis Model (k = 234; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 1.9496 1.3963 21 no study_ID no
## sigma^2.2 0.0001 0.0074 78 no spp_name_phylo yes
## sigma^2.3 0.0925 0.3042 234 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 232) = 1579.6588, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 232) = 0.1117, p-val = 0.7385
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.2211 0.3139 -0.7043 232 0.4819 -0.8397 0.3974
## SSD_index -0.2015 0.6028 -0.3342 232 0.7385 -1.3890 0.9861
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # lnCVR
MLMR_mods_pers_SSD_bird_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"),
test = "t", data = pers_new_bird)## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_bird_lncvr##
## Multivariate Meta-Analysis Model (k = 234; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0000 0.0003 21 no study_ID no
## sigma^2.2 0.0029 0.0537 78 no spp_name_phylo yes
## sigma^2.3 0.0000 0.0003 234 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 232) = 244.9667, p-val = 0.2670
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 232) = 0.6126, p-val = 0.4346
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0428 0.0361 1.1851 232 0.2372 -0.0283 0.1139
## SSD_index 0.1023 0.1307 0.7827 232 0.4346 -0.1553 0.3599
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.5.4 Fish
SSD for aggression and boldness.
# 3. FISH
# subset by trait type
# a. aggression
pers_new_fish_aggression <- as.data.frame(pers_new %>%
filter(personality_trait == "aggression") %>%
filter(taxo_group == "fish"))
# b. boldness
pers_new_fish_bold <- as.data.frame(pers_new %>%
filter(personality_trait == "boldness") %>%
filter(taxo_group == "fish"))
# phylo
phylo_vcv_fish <- phylo_vcv[[2]]Aggression:
# a. aggression
# SMD
MLMR_mods_pers_SSD_fish_aggression_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"),
test = "t", data = pers_new_fish_aggression)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.## Warning: Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_fish_aggression_SMD##
## Multivariate Meta-Analysis Model (k = 93; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0194 0.1395 16 no study_ID no
## sigma^2.2 0.3329 0.5770 13 no spp_name_phylo yes
## sigma^2.3 0.1704 0.4128 93 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 91) = 334.1728, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 91) = 0.2301, p-val = 0.6326
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.1643 0.3987 -0.4120 91 0.6813 -0.9562 0.6277
## SSD_index 0.2659 0.5544 0.4797 91 0.6326 -0.8352 1.3671
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # lnCVR
MLMR_mods_pers_SSD_fish_aggression_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"),
test = "t", data = pers_new_fish_aggression)## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_fish_aggression_lncvr##
## Multivariate Meta-Analysis Model (k = 93; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0210 0.1450 16 no study_ID no
## sigma^2.2 0.0000 0.0021 13 no spp_name_phylo yes
## sigma^2.3 0.0000 0.0012 93 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 91) = 68.2701, p-val = 0.9640
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 91) = 0.1495, p-val = 0.6999
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.1163 0.0597 -1.9483 91 0.0545 -0.2348 0.0023 .
## SSD_index -0.1323 0.3423 -0.3866 91 0.6999 -0.8122 0.5476
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Boldness:
# SMD
MLMR_mods_pers_SSD_fish_bold_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"),
test = "t", data = pers_new_fish_bold)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.## Warning: Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_fish_bold_SMD##
## Multivariate Meta-Analysis Model (k = 171; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1717 0.4143 23 no study_ID no
## sigma^2.2 0.0279 0.1671 12 no spp_name_phylo yes
## sigma^2.3 0.1634 0.4042 171 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 169) = 614.1157, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 169) = 0.3196, p-val = 0.5726
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.1883 0.1764 1.0678 169 0.2871 -0.1599 0.5366
## SSD_index -0.2571 0.4548 -0.5653 169 0.5726 -1.1550 0.6408
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # lnCVR
MLMR_mods_pers_SSD_fish_bold_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"),
test = "t", data = pers_new_fish_bold) ## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_fish_bold_lncvr##
## Multivariate Meta-Analysis Model (k = 171; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0445 0.2109 23 no study_ID no
## sigma^2.2 0.0184 0.1356 12 no spp_name_phylo yes
## sigma^2.3 0.0881 0.2968 171 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 169) = 395.3375, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 169) = 0.3068, p-val = 0.5804
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0063 0.1310 0.0482 169 0.9616 -0.2523 0.2649
## SSD_index -0.1507 0.2721 -0.5539 169 0.5804 -0.6878 0.3864
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.5.5 Inverts
SSD for activity and boldness.
# 4. INVERTS
# subset dataset
# a. activity
invert_activity <- as.data.frame(pers_new %>%
filter(personality_trait == "activity") %>%
filter(taxo_group == "invertebrate"))
# b. boldness
invert_bold <- as.data.frame(pers_new %>%
filter(personality_trait == "boldness") %>%
filter(taxo_group == "invertebrate"))
# phylo
phylo_vcv_invert <- phylo_vcv[[3]]Activity:
# rerun models
# a. activity
# SMD
MLMR_mods_pers_SSD_invert_activity_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"),
test = "t", data = invert_activity)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.## Warning: Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_invert_activity_SMD##
## Multivariate Meta-Analysis Model (k = 165; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 2.1874 1.4790 18 no study_ID no
## sigma^2.2 0.0001 0.0104 16 no spp_name_phylo yes
## sigma^2.3 0.1562 0.3952 165 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 163) = 1081.7241, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 163) = 0.7140, p-val = 0.3993
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.3479 0.3670 0.9480 163 0.3445 -0.3767 1.0725
## SSD_index -0.6862 0.8120 -0.8450 163 0.3993 -2.2896 0.9173
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # lnCVR
MLMR_mods_pers_SSD_invert_activity_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"),
test = "t", data = invert_activity) ## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data.
## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## Rows with NAs omitted from model fitting. MLMR_mods_pers_SSD_invert_activity_lncvr##
## Multivariate Meta-Analysis Model (k = 165; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1193 0.3454 18 no study_ID no
## sigma^2.2 0.0000 0.0035 16 no spp_name_phylo yes
## sigma^2.3 0.0591 0.2431 165 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 163) = 486.7410, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 163) = 0.4392, p-val = 0.5085
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0278 0.1074 -0.2589 163 0.7960 -0.2400 0.1843
## SSD_index 0.2685 0.4052 0.6627 163 0.5085 -0.5315 1.0685
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Boldness:
# SMD
MLMR_mods_pers_SSD_invert_bold_SMD <- rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"),
test = "t", data = invert_bold)## Warning in rma.mv(SMD_yi_flip ~ SSD_index, V = SMD_vi, random = list(~1 | :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data. MLMR_mods_pers_SSD_invert_bold_SMD##
## Multivariate Meta-Analysis Model (k = 164; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0822 0.2867 23 no study_ID no
## sigma^2.2 0.0000 0.0020 23 no spp_name_phylo yes
## sigma^2.3 0.0650 0.2550 164 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 162) = 513.4222, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 162) = 0.1533, p-val = 0.6959
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0985 0.0823 1.1967 162 0.2332 -0.0640 0.2611
## SSD_index 0.1313 0.3354 0.3915 162 0.6959 -0.5310 0.7936
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 # lnCVR
MLMR_mods_pers_SSD_invert_bold_lncvr <- rma.mv(CVR_yi ~ SSD_index, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_invert), control=list(optimizer="optim"),
test = "t", data = invert_bold) ## Warning in rma.mv(CVR_yi ~ SSD_index, V = CVR_vi, random = list(~1 | study_ID, :
## There are rows/columns in the 'R' matrix for 'spp_name_phylo' for which there
## are no data. MLMR_mods_pers_SSD_invert_bold_lncvr##
## Multivariate Meta-Analysis Model (k = 164; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0383 0.1957 23 no study_ID no
## sigma^2.2 0.0000 0.0015 23 no spp_name_phylo yes
## sigma^2.3 0.0378 0.1945 164 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 162) = 380.7049, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 162) = 0.0013, p-val = 0.9716
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0145 0.0607 -0.2384 162 0.8119 -0.1344 0.1055
## SSD_index -0.0089 0.2503 -0.0356 162 0.9716 -0.5032 0.4854
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 18.6 Multiple testing
Now we can extract the p-values from our intercept models, personality trait models, and SSD subset models to adjust p-values using the false discovery rate method. This method uses the p.adjust function to adjust p-values to account for multiple testing.
## Extract p-values from SSD subset models (reported in main text)
# list
p.SMD_SSD <- c(0.69, 0.03, 0.86, 0.04, 0.34, 0.25, 0.99, 0.85, 0.48, 0.74, 0.68, 0.63, 0.29, 0.57, 0.34, 0.40, 0.23, 0.70)
p.lnCVR_SSD <- c(0.60, 0.72, 0.52, 0.91, 0.82, 0.15, 0.69, 0.61, 0.24, 0.43, 0.05, 0.70, 0.96, 0.58, 0.80, 0.51, 0.81, 0.97)
# p adjustment on our 3 hypothesis-testing models
#SMD
p.adjust(p = c(p.SMD_intercept, p.SMD_pers, p.SMD_SSD), method = "fdr") ## bird fish invertebrate mammal reptilia
## 0.8533333 0.8533333 0.8533333 0.8708500 0.8533333
## bird1 bird2 bird3 bird4 bird5
## 0.8708500 0.8533333 0.8533333 0.8708500 0.4800000
## fish1 fish2 fish3 fish4 fish5
## 0.8708500 0.8533333 0.8708500 0.8533333 0.8708500
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5
## 0.8533333 0.8533333 0.8533333 0.9195744 0.8533333
## mammal1 mammal2 mammal3 mammal4 mammal5
## 0.8533333 0.8708500 0.8533333 0.9292872 0.8708500
## reptilia1 reptilia2 reptilia3 reptilia4 reptilia5
## 0.9195744 0.8533333 0.8533333 0.4800000 0.9900000
##
## 0.8708500 0.4800000 0.9195744 0.4800000 0.8533333
##
## 0.8533333 0.9900000 0.9195744 0.8533333 0.8708500
##
## 0.8708500 0.8708500 0.8533333 0.8708500 0.8533333
##
## 0.8533333 0.8533333 0.8708500 #lncvr
p.adjust(p = c(p.lnCVR_intercept, p.lnCVR_pers, p.lnCVR_SSD), method = "fdr") ## bird fish invertebrate mammal reptilia
## 0.9513857 0.9513857 0.9513857 0.9513857 0.9513857
## bird1 bird2 bird3 bird4 bird5
## 0.9513857 0.9513857 0.9513857 0.4948879 0.9513857
## fish1 fish2 fish3 fish4 fish5
## 0.9513857 0.8505594 0.9513857 0.9513857 0.4948879
## invertebrate1 invertebrate2 invertebrate3 invertebrate4 invertebrate5
## 0.9513857 0.9513857 0.9513857 0.9513857 0.9513857
## mammal1 mammal2 mammal3 mammal4 mammal5
## 0.9513857 0.9513857 0.9513857 0.9513857 0.9513857
## reptilia1 reptilia2 reptilia3 reptilia4 reptilia5
## 0.9513857 0.7414563 0.9513857 0.9513857 0.9513857
##
## 0.9513857 0.9513857 0.9513857 0.9513857 0.9513857
##
## 0.9513857 0.9513857 0.9513857 0.9513857 0.9513857
##
## 0.7414563 0.9513857 0.9700000 0.9513857 0.9513857
##
## 0.9513857 0.9513857 0.9700000 # these p-values are in the order presented in tables, so easy to replace old p-values with new ones9 Exploratory analyses
We collected some additional information from the literature (mating system) and from studies that we expected would influence sex differences. These analyses are strictly exploratory and just compare categorical moderator terms.
9.1 mating system
Do effect sizes from monogamous or multiply-mating species differ? Model summaries presented in Supplementary Table S3.
# what have we got to work with?
pers_new %>%
group_by(taxo_group, mating_system) %>%
filter(!is.na(mating_system))%>%
summarise(n = n(), studies = length(unique(study_ID)), species = length(unique(spp_name_phylo))) # make a table of numbers## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.# reload model output
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMR_models_pers_mating_system <- meta_model_fits(pers_new, phylo_vcv, type = "pers_mate")
saveRDS(MLMR_models_pers_mating_system, "./output/MLMR_models_pers_mating_system")
} else{
MLMR_models_pers_mating_system <- readRDS("./output/MLMR_models_pers_mating_system")
}
# Extract the SMD and lnCVR results
smd_mods_mating_system <- MLMR_models_pers_mating_system["SMD",]
lnCVR_mods_mating_system <- MLMR_models_pers_mating_system["lnCVR",]9.2 age
Do effect sizes from adults (sexually mature) or juveniles differ? Model summaries presented in Supplementary Table S4
# make a table
data.frame(pers_new %>%
group_by(taxo_group, age) %>%
summarise(n= n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.# reload model output:
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMR_models_pers_age <- meta_model_fits(pers_new, phylo_vcv, type = "age")
saveRDS(MLMR_models_pers_age, "./output/MLMR_models_pers_age")
} else{
MLMR_models_pers_age <- readRDS("./output/MLMR_models_pers_age")
}
# Extract the SMD and lnCVR results
smd_mods_pers_age <- MLMR_models_pers_age["SMD",]
lnCVR_mods_pers_age <- MLMR_models_pers_age["lnCVR",]9.3 population
Do effect sizes from wild animals or lab animals differ? Model summaries presented in Supplementary Table S5.
# table
data.frame(pers_new %>%
group_by(taxo_group, population) %>%
summarise(n = n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.# reload model output:
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMR_models_pers_pop <- meta_model_fits(pers_new, phylo_vcv, type = "pop")
saveRDS(MLMR_models_pers_pop, "./output/MLMR_models_pers_pop")
} else{
MLMR_models_pers_pop <- readRDS("./output/MLMR_models_pers_pop")
}
# Extract the SMD and lnCVR results
smd_mods_pers_pop <- MLMR_models_pers_pop["SMD",]
lnCVR_mods_pers_pop <- MLMR_models_pers_pop["lnCVR",]9.4 study environment
Do effect sizes collected in the wild or the lab differ? Model summaries presented in Supplementary Table S6.
# table
data.frame(pers_new %>%
group_by(taxo_group, study_environment) %>%
summarise(n = n(), N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument.# reload model output:
rerun_models == FALSE## [1] TRUE if(rerun_models == TRUE){
MLMR_models_pers_environ <- meta_model_fits(pers_new, phylo_vcv, type = "environ")
saveRDS(MLMR_models_pers_environ, "./output/MLMR_models_pers_environ")
} else{
MLMR_models_pers_environ <- readRDS("./output/MLMR_models_pers_environ")
}
# Extract the SMD and lnCVR results
smd_mods_pers_enviro <- MLMR_models_pers_environ["SMD",]
lnCVR_mods_pers_enviro <- MLMR_models_pers_environ["lnCVR",]9.5 study type
Do effect sizes from observational or experimental study design differ? Model summaries presented in Supplementary Table S7.
# Let's see what we have to work with
data.frame(pers_new %>%
group_by(taxo_group, study_type) %>%
summarise(N_spp = length(unique(spp_name_phylo)), N_studies = length(unique(study_ID))))## `summarise()` has grouped output by 'taxo_group'. You can override using the `.groups` argument. # inverts only have experimental observations, so need to exclude inverts from this analysis
# because our phylo_vcv matrix is in a list that is hard to drop elements from, let's just run each model individually
# 1. Mammals
# Subset data
pers_new_mammal <- as.data.frame(pers_new %>%
filter(taxo_group == "mammal"))
# Run models with just study type as moderator:
#SMD
MLMR_mods_pers_studytype_mammal_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal)
MLMR_mods_pers_studytype_mammal_SMD##
## Multivariate Meta-Analysis Model (k = 674; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.1051 0.3242 61 no study_ID no
## sigma^2.2 0.0094 0.0971 45 no spp_name_phylo yes
## sigma^2.3 0.1570 0.3963 674 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 672) = 2198.5222, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 672) = 9.6851, p-val = 0.0019
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0067 0.0913 -0.0734 672 0.9415 -0.1860 0.1726
## study_typeobservation 0.4092 0.1315 3.1121 672 0.0019 0.1510 0.6674
##
## intrcpt
## study_typeobservation **
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_studytype_mammal_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_mammal), control=list(optimizer="optim"),
test = "t", data = pers_new_mammal)
MLMR_mods_pers_studytype_mammal_lncvr##
## Multivariate Meta-Analysis Model (k = 674; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0356 0.1888 61 no study_ID no
## sigma^2.2 0.0436 0.2088 45 no spp_name_phylo yes
## sigma^2.3 0.0338 0.1837 674 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 672) = 1061.6294, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 672) = 0.5661, p-val = 0.4521
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0298 0.1361 0.2188 672 0.8269 -0.2375 0.2970
## study_typeobservation 0.0771 0.1024 0.7524 672 0.4521 -0.1240 0.2781
##
## intrcpt
## study_typeobservation
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 2. BIRDS
# subset dataset
pers_new_bird <- as.data.frame(pers_new %>%
filter(taxo_group == "bird"))
# rerun models
#SMD
MLMR_mods_pers_studytype_bird_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"),
test = "t", data = pers_new_bird)
MLMR_mods_pers_studytype_bird_SMD##
## Multivariate Meta-Analysis Model (k = 480; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.6650 0.8155 50 no study_ID no
## sigma^2.2 0.0000 0.0038 106 no spp_name_phylo yes
## sigma^2.3 0.1203 0.3469 480 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 478) = 2378.8387, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 478) = 1.0068, p-val = 0.3162
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0210 0.1515 -0.1389 478 0.8896 -0.3187 0.2766
## study_typeobservation -0.2540 0.2532 -1.0034 478 0.3162 -0.7515 0.2434
##
## intrcpt
## study_typeobservation
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_studytype_bird_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_bird), control=list(optimizer="optim"),
test = "t", data = pers_new_bird)
MLMR_mods_pers_studytype_bird_lncvr##
## Multivariate Meta-Analysis Model (k = 480; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.2537 0.5036 50 no study_ID no
## sigma^2.2 0.0001 0.0086 106 no spp_name_phylo yes
## sigma^2.3 0.3766 0.6137 480 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 478) = 3186.9869, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 478) = 1.7076, p-val = 0.1919
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0436 0.1109 0.3929 478 0.6946 -0.1744 0.2615
## study_typeobservation -0.2451 0.1876 -1.3067 478 0.1919 -0.6138 0.1235
##
## intrcpt
## study_typeobservation
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 3. FISH
# subset dataset
pers_new_fish <- as.data.frame(pers_new %>%
filter(taxo_group == "fish"))
# rerun models
#SMD
MLMR_mods_pers_studytype_fish_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"),
test = "t", data = pers_new_fish)
MLMR_mods_pers_studytype_fish_SMD##
## Multivariate Meta-Analysis Model (k = 490; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.5997 0.7744 44 no study_ID no
## sigma^2.2 0.0440 0.2098 22 no spp_name_phylo yes
## sigma^2.3 0.1091 0.3303 490 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 488) = 1523.5807, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 488) = 0.0188, p-val = 0.8911
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.1825 0.2106 0.8669 488 0.3864 -0.2312 0.5962
## study_typeobservation -0.0657 0.4798 -0.1369 488 0.8911 -1.0085 0.8771
##
## intrcpt
## study_typeobservation
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_studytype_fish_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_fish), control=list(optimizer="optim"),
test = "t", data = pers_new_fish)
MLMR_mods_pers_studytype_fish_lncvr##
## Multivariate Meta-Analysis Model (k = 490; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0352 0.1875 44 no study_ID no
## sigma^2.2 0.0017 0.0408 22 no spp_name_phylo yes
## sigma^2.3 0.1094 0.3307 490 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 488) = 1122.6615, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 488) = 0.0722, p-val = 0.7882
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt -0.0025 0.0536 -0.0460 488 0.9633 -0.1078 0.1029
## study_typeobservation -0.0399 0.1486 -0.2687 488 0.7882 -0.3319 0.2520
##
## intrcpt
## study_typeobservation
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# 4. Reptiles
# subset dataset
pers_new_reptile <- as.data.frame(pers_new %>%
filter(taxo_group == "reptilia"))
# phylo
phylo_vcv_reptile <- phylo_vcv[[5]]
# rerun models
#SMD
MLMR_mods_pers_studytype_rep_SMD <- rma.mv(SMD_yi_flip ~ study_type, V = SMD_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_reptile), control=list(optimizer="optim"),
test = "t", data = pers_new_reptile)
MLMR_mods_pers_studytype_rep_SMD ##
## Multivariate Meta-Analysis Model (k = 95; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0000 0.0008 11 no study_ID no
## sigma^2.2 0.0730 0.2702 10 no spp_name_phylo yes
## sigma^2.3 0.0426 0.2063 95 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 93) = 159.2776, p-val < .0001
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 93) = 3.4462, p-val = 0.0666
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.1334 0.1531 0.8715 93 0.3857 -0.1706 0.4374
## study_typeobservation -0.5085 0.2739 -1.8564 93 0.0666 -1.0524 0.0354
##
## intrcpt
## study_typeobservation .
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 #lnCVR
MLMR_mods_pers_studytype_rep_lncvr <- rma.mv(CVR_yi ~ study_type, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylo_vcv_reptile), control=list(optimizer="optim"),
test = "t", data = pers_new_reptile)
MLMR_mods_pers_studytype_rep_lncvr##
## Multivariate Meta-Analysis Model (k = 95; method: REML)
##
## Variance Components:
##
## estim sqrt nlvls fixed factor R
## sigma^2.1 0.0000 0.0003 11 no study_ID no
## sigma^2.2 0.0000 0.0012 10 no spp_name_phylo yes
## sigma^2.3 0.0000 0.0003 95 no obs no
##
## Test for Residual Heterogeneity:
## QE(df = 93) = 58.0505, p-val = 0.9983
##
## Test of Moderators (coefficient 2):
## F(df1 = 1, df2 = 93) = 0.1931, p-val = 0.6613
##
## Model Results:
##
## estimate se tval df pval ci.lb ci.ub
## intrcpt 0.0427 0.0420 1.0173 93 0.3116 -0.0407 0.1261
## study_typeobservation -0.0621 0.1413 -0.4395 93 0.6613 -0.3426 0.2184
##
## intrcpt
## study_typeobservation
##
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 110 Sensitivity analyses - Dependency matrix models
We need to refit our 3 main models accounting for any dependency resulting from the same traits measured on the same animals (likely a big source of non-independence) and any other shared covariance. We added the D matrices to the residual variance matrix as opposed to the sampling covariance. We chose to set 3 different levels of dependency (rho): 0.3, 0.5 an 0.8.
Model summaries are also presented in Supplementary Tables S8-S13.
# Create the dependency matrices; try 3 levels of rho = 0.3, 0.5, 0.8
pers_new <- data.frame(pers_new %>%
group_by(taxo_group) %>%
mutate(depend_n = paste0(study_ID, "_", depend)))
split_taxa <- split(pers_new, pers_new$taxo_group)
# 0.3 rho:
D_matrices_0.3 <- lapply(split_taxa, function(x) make_VCV_matrix(x, V = x$SMD_vi, cluster = "depend_n",
obs = "obs", type = "cor", rho = 0.3))
# 0.5 rho:
D_matrices_0.5 <- lapply(split_taxa, function(x) make_VCV_matrix(x, V = x$SMD_vi, cluster = "depend_n",
obs = "obs", type = "cor", rho = 0.5))
# 0.8 rho:
D_matrices_0.8 <- lapply(split_taxa, function(x) make_VCV_matrix(x, V = x$SMD_vi, cluster = "depend_n",
obs = "obs", type = "cor", rho = 0.8))10.1 Intercept-only models with D matrices
Model output is presented in Supplementary Tables S8-S10 in the Supporting Information.
# 1. Intercept only models
# rho = 0.3
int_0.3 <- fit_int_MLMAmodel_D(pers_new, phylo_vcv, D_matrices_0.3)
smd_mods_D_0.3 <- int_0.3[["SMD"]]
lnCVR_mods_D_0.3 <- int_0.3[["lnCVR"]]
# prediction intervals
MLMA_estimates_SMD_D_0.3 <- plyr::ldply(lapply(smd_mods_D_0.3,
function(x) print(mod_results(x, mod = "Int"))))
MLMA_estimates_lnCVR_D_0.3 <- plyr::ldply(lapply(lnCVR_mods_D_0.3,
function(x) print(mod_results(x, mod = "Int"))))
# rho = 0.5
int_0.5 <- fit_int_MLMAmodel_D(pers_new, phylo_vcv, D_matrices_0.5)
smd_mods_D_0.5 <- int_0.5[["SMD"]]
lnCVR_mods_D_0.5 <- int_0.5[["lnCVR"]]
# prediction intervals
MLMA_estimates_SMD_D_0.5 <- plyr::ldply(lapply(smd_mods_D_0.5,
function(x) print(mod_results(x, mod = "Int"))))
MLMA_estimates_lnCVR_D_0.5 <- plyr::ldply(lapply(lnCVR_mods_D_0.5,
function(x) print(mod_results(x, mod = "Int"))))
# rho = 0.8
int_0.8 <- fit_int_MLMAmodel_D(pers_new, phylo_vcv, D_matrices_0.8)
smd_mods_D_0.8 <- int_0.8[["SMD"]]
lnCVR_mods_D_0.8 <- int_0.8[["lnCVR"]]
# prediction intervals
MLMA_estimates_SMD_D_0.8 <- plyr::ldply(lapply(smd_mods_D_0.8,
function(x) print(mod_results(x, mod = "Int"))))
MLMA_estimates_lnCVR_D_0.8 <- plyr::ldply(lapply(lnCVR_mods_D_0.8,
function(x) print(mod_results(x, mod = "Int"))))10.2 Personality trait models with D matrices
Model output is presented in Supplementary Tables S11-S13 in the Supporting Information.
# 2. Personality Trait models
# rho = 0.3
pers_0.3 <- fit_int_MLMAmodel_D_pers(pers_new, phylo_vcv, D_matrices_0.3)
smd_mods_D_pers_0.3 <- pers_0.3[["SMD"]]
lnCVR_mods_D_pers_0.3 <- pers_0.3[["lnCVR"]]
# prediction intervals
MLMA_estimates_SMD_pers_D_0.3 <- plyr::ldply(lapply(smd_mods_D_pers_0.3,
function(x) print(mod_results(x, mod = "personality_trait"))))
MLMA_estimates_lnCVR_pers_D_0.3 <- plyr::ldply(lapply(lnCVR_mods_D_pers_0.3,
function(x) print(mod_results(x, mod = "personality_trait"))))
# rho = 0.5
pers_0.5 <- fit_int_MLMAmodel_D_pers(pers_new, phylo_vcv, D_matrices_0.5)
smd_mods_D_pers_0.5 <- pers_0.5[["SMD"]]
lnCVR_mods_D_pers_0.5 <- pers_0.5[["lnCVR"]]
# prediction intervals
MLMA_estimates_SMD_pers_D_0.5 <- plyr::ldply(lapply(smd_mods_D_pers_0.5,
function(x) print(mod_results(x, mod = "personality_trait"))))
MLMA_estimates_lnCVR_pers_D_0.5 <- plyr::ldply(lapply(lnCVR_mods_D_pers_0.5,
function(x) print(mod_results(x, mod = "personality_trait"))))
# rho = 0.8
pers_0.8 <- fit_int_MLMAmodel_D_pers(pers_new, phylo_vcv, D_matrices_0.8)
smd_mods_D_pers_0.8 <- pers_0.8[["SMD"]]
lnCVR_mods_D_pers_0.8 <- pers_0.8[["lnCVR"]]
# prediction intervals
MLMA_estimates_SMD_pers_D_0.8 <- plyr::ldply(lapply(smd_mods_D_pers_0.8,
function(x) print(mod_results(x, mod = "personality_trait"))))
MLMA_estimates_lnCVR_pers_D_0.8 <- plyr::ldply(lapply(lnCVR_mods_D_pers_0.8,
function(x) print(mod_results(x, mod = "personality_trait")))) 10.3 Personality * SSD models
These models were just to check since we don't really interpret the interaction models.
# 3. Pers Trait * SSD models
# just use the full interaction models here since this is just a check
# won't bother with prediction intervals here since these models aren't really for interpretation
# rho = 0.3
ssd_0.3 <- fit_int_MLMAmodel_D_pers_ssd(pers_new, phylo_vcv, D_matrices_0.3)
split_taxa <- split(pers_new, pers_new$taxo_group)
smd_mods_D_pers_ssd_0.3 <- ssd_0.3[["SMD"]]
lnCVR_mods_D_pers_ssd_0.3 <- ssd_0.3[["lnCVR"]]
# rho = 0.5
ssd_0.5 <- fit_int_MLMAmodel_D_pers_ssd(pers_new, phylo_vcv, D_matrices_0.5)
split_taxa <- split(pers_new, pers_new$taxo_group)
smd_mods_D_pers_ssd_0.5 <- ssd_0.5[["SMD"]]
lnCVR_mods_D_pers_ssd_0.5 <- ssd_0.5[["lnCVR"]]
# rho = 0.8
ssd_0.8 <- fit_int_MLMAmodel_D_pers_ssd(pers_new, phylo_vcv, D_matrices_0.8)
split_taxa <- split(pers_new, pers_new$taxo_group)
smd_mods_D_pers_ssd_0.8 <- ssd_0.8[["SMD"]]
lnCVR_mods_D_pers_ssd_0.8 <- ssd_0.8[["lnCVR"]] 11 Sensitivity analyses - Publication Bias
We can use: 1) funnel plots to look for asymmetry across all effect sizes for both SMD and lnCVR, and 2) Egger's test which performs a regression test on our funnel plots ... but is not useful when there is high heterogeneity NOT caused by publication bias (which is the case for our data).
Since our data has very high heterogeneity, we instead included the inverse of the 'effective sample size' as a moderator term in our full model (personality trait model) to see if study precision is driving effect size patterns. The logic here is that studies with low or high precision can have a significant influence and so including precision as a moderator will allow us to see if precision is significant (and which direction). See Nakagawa et al. 2021 for more info (reference in main text).
Model summaries are presented in Supplementary Table S18.
### NEW METHOD OF PUBLICATION BIAS FROM NAKAGAWA ET AL 2021 - PREPRINT
# calculating the inverse of the 'effective sample size' to account for unbalanced sampling
pers_new$inv_n_tilda <- with(pers_new, ((female_n + male_n)/(female_n*male_n)))
pers_new$sqrt_inv_n_tilda <- with(pers_new, (sqrt(inv_n_tilda))) # use this in the model
if(rerun_models == TRUE){
MLMR_models_pers_pubbias <- meta_model_fits(pers_new, phylo_vcv, type = "pubbias")
saveRDS(MLMR_models_pers_pubbias, "./output/MLMR_models_pers_pubbias")
} else{
MLMR_models_pers_pubbias <- readRDS("./output/MLMR_models_pers_pubbias")
}
# Extract the SMD and lnCVR results
smd_mods_pubbias <- MLMR_models_pers_pubbias["SMD",]
lnCVR_mods_pubbias <- MLMR_models_pers_pubbias["lnCVR",] 12 Exploratory analyses - Heterogamety and taxo group
There was a trend for male mammals to be more variable than females and female birds to be more variable than males, for some of the five personality traits. To better compare the direction of these effect sizes we decided post hoc to conduct an exploratory analysis with personality trait type and taxonomic group as moderator terms to compare birds and mammals (males homogametic or heterogametic, respectively). To do this, we first combined the bird and mammal phylo correlation matrices together (assuming no phylo heritability across the taxo groups - since phylo did not really explain heterogeneity it shouldn't matter). We then created an interaction MLMR model with personality trait * taxa (no intercept) to get slope estimates for each of the traits for mammals and birds seperately.
From this model, we then compared each of the five traits for birds and mammals using a post hoc Tukey pairwise comparison to test whether birds and mammals were significantly different from each other.
Model summaries are presented in Supplementary Table S19.
# install packages to make diagonal matrix and to make multiple comparisons
library(multcomp)
library(Matrix)
# Create block diag phylogeny
phylogeny <- Matrix::bdiag(phylo_vcv_bird, phylo_vcv_mammal) # use this as the phylo vcv in the model
# needs to have colnames for use in random effects model
dimnames(phylogeny) <- Map(c, dimnames(phylo_vcv_bird), dimnames(phylo_vcv_mammal))
# only include bird and mammal data
pers_new_contrast <- as.data.frame(pers_new %>%
filter(taxo_group =="mammal" | taxo_group == "bird") %>%
mutate(sp_pers = interaction(personality_trait,taxo_group)))
# 1. intercept only model
contrast_birdmammal_lncvr_int <- rma.mv(CVR_yi ~ taxo_group, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylogeny), control=list(optimizer="optim"),
test = "t", data = pers_new_contrast)
# 2. personality trait model
# creating the model - with pers trait and taxo group as mods
#lnCVR model only
contrast_birdmammal_lncvr <- rma.mv(CVR_yi ~ personality_trait*taxo_group, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylogeny), control=list(optimizer="optim"),
test = "t", data = pers_new_contrast)
# model with interaction only to check output of model above
contrast_birdmammal_lncvr_2 <- rma.mv(CVR_yi ~ sp_pers -1, V = CVR_vi,
random = list(~1|study_ID, ~1|spp_name_phylo, ~1|obs),
R = list(spp_name_phylo=phylogeny), control=list(optimizer="optim"),
test = "t", data = pers_new_contrast)
# multiple comparison using Tukey test
summary(glht(contrast_birdmammal_lncvr, linfct = cbind(contrMat(rep(1:10), type = "Tukey"))), test=adjusted("fdr"))
# here we are only interested in the comparisons between mammals and birds, so: 1-6 (activity), 2-7 (aggression), 3-8 (boldness), 4-9 (exploration), and 5-10 (sociality)13 Plots
13.1 Orchard plots of effect sizes from personality trait models
These plots use the orchaRd package to generate pretty plots where each effect size (k) is a point on the plot.
13.1.1 lnCVR
# create objects of each of the models first
# lnCVR
# Bird lnCVR
bird_lncvr <- orchard_plot(lnCVR_mods_pers[[1]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
# Fish lnCVR
fish_lncvr <- orchard_plot(lnCVR_mods_pers[[2]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
# Invert lnCVR
invert_lncvr<- orchard_plot(lnCVR_mods_pers[[3]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
# Mammal lnCVR
mammal_lncvr <- orchard_plot(lnCVR_mods_pers[[4]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")
# Reptile lnCVR
reptile_lncvr <- orchard_plot(lnCVR_mods_pers[[5]], mod = "personality_trait", xlab = "log Coefficient of Variance (lnCVR)", angle = 45, alpha = 0.5, transfm = "none")13.1.2 SMD
# Bird SMD
bird_SMD <- orchard_plot(smd_mods_pers[[1]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
# Fish SMD
fish_SMD <- orchard_plot(smd_mods_pers[[2]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
# Invert SMD
invert_SMD<- orchard_plot(smd_mods_pers[[3]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
# Mammal SMD
mammal_SMD <- orchard_plot(smd_mods_pers[[4]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")
# Reptile SMD
reptile_SMD <- orchard_plot(smd_mods_pers[[5]], mod = "personality_trait", xlab = "Standardised mean difference", angle = 45, alpha = 0.5, transfm = "none")Putting the SMD and lnCVR plots together
Endotherms:
# window size for orchard plots
# the precision guides on the plots are a bit ugly, collect them to the side and crop them out
## Mammals
# dev.new(width=8,height=7,noRStudioGD = TRUE)
# place plots together
mammal_SMD / mammal_lncvr 
# ggsave("./figs/finished figs/mammal_effects.tiff", width = 8, height = 7, units = "in") #save image
## Birds
bird_SMD / bird_lncvr
# ggsave("./figs/finished figs/bird_effects.tiff", width = 8, height = 7, units = "in") #save imageEctotherms:
# window size a bit smaller for these guys
# dev.new(width=8,height=5,noRStudioGD = TRUE)
## Reptiles and amphibians
reptile_SMD / reptile_lncvr / plot_layout(guides = 'collect')
# ggsave("~/Documents/GitHub/sex_meta/figs/finished figs/rep_effects.tiff", width = 8, height = 5, units = "in")
## Fish
fish_SMD / fish_lncvr / plot_layout(guides = 'collect')
# ggsave("~/Documents/GitHub/sex_meta/figs/finished figs/fish_effects.tiff", width = 8, height = 5, units = "in")
## Invertebrates
invert_SMD / invert_lncvr / plot_layout(guides = 'collect')
# ggsave("~/Documents/GitHub/sex_meta/figs/finished figs/invert_effects.tiff", width = 8, height = 5, units = "in")The precision guides will get cropped out when joining the orchard plots and phylogenies together.
13.2 Phylogenetic trees with heatmaps
Using ggtree to plot lots of complex data onto phylogenetic trees see: https://guangchuangyu.github.io/ggtree-book/chapter-ggtree.html for more information about using ggtree
# install ggtree using this method:
source("https://bioconductor.org/biocLite.R")
BiocManager::install("ggtree")##
## The downloaded binary packages are in
## /var/folders/0b/pxghylq157gfhs1vrzdpx2gc0000gq/T//RtmpXG0rME/downloaded_packages library(ggtree)
# load organised SSD data using figs_data.csv
figs_data <- read.csv("./data/figs_data.csv", stringsAsFactors = FALSE)13.3 bird tree
# subset dataset to include only birds
bird_data <- as.data.frame(figs_data %>%
filter(taxo_group == "bird"))
# setting up the basic tree structure
# load tree
birdtree <- read.tree("./trees/bird_species.nwk")
# prune tree to get rid of species we no longer have data for
pruned.birdtree <- drop.tip(birdtree, setdiff(birdtree$tip.label, bird_data$spp_name_phylo))
# remove underscores from tip labels
pruned.birdtree$tip.label = gsub("_", " ", pruned.birdtree$tip.label)
# remove underscores from species name in our species data list
bird_data$spp_name_phylo = gsub("_", " ", bird_data$spp_name_phylo)
# set row names
row.names(bird_data) <- bird_data$spp_name_phylo
# define objects for the plot
species <- pruned.birdtree$tip.label
rownames(bird_data) <- pruned.birdtree$tip.label
# set window size
# dev.new(width=8, height=8,noRStudioGD = TRUE) #opens quartz window of set size
# now need to make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
# subset dataset
pers_bird <- as.data.frame(pers_new %>%
filter(taxo_group == "bird"))
# make this a matrix-style dataframe
pers_bird <- data.frame(pers_bird %>%
group_by(spp_name_phylo, personality_trait) %>%
summarise(n = n()))## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument. # remove underscores species names
pers_bird$spp_name_phylo = gsub("_", " ", pers_bird$spp_name_phylo)
# create matrix
pers_bird <- data.frame(pers_bird %>%
spread(personality_trait, n, fill = 0))
# set species name as row name for matrix
row.names(pers_bird) <- pers_bird$spp_name_phylo
pers_bird <- pers_bird[,2:6]
# matrix
birds_matrix <- data.matrix(pers_bird)
# FINAL TREE
# making the tree
p_b1 <- ggtree(pruned.birdtree, size = 0.3, layout = 'circular', branch.length = 'none') %<+% bird_data +
xlim(-40, NA) +
geom_tippoint(aes(color=SSD_index)) +
scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") +
geom_tiplab2(size = 2.2, offset = 4, colour = "black", fontface = "italic") +
theme(legend.position = 'right')## Warning: `data_frame()` was deprecated in tibble 1.1.0.
## Please use `tibble()` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.## Warning: `mutate_()` was deprecated in dplyr 0.7.0.
## Please use `mutate()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated. # adding heatmap of traits
p_b2 <- gheatmap(p_b1, birds_matrix, offset=68, width=2, low = "white", high = "mediumseagreen", color=NULL,
colnames=T, colnames_angle = 60, colnames_offset_y = .1, colnames_offset_x = .2) +
theme(plot.tag = element_text(size = 2, face = "bold"),
legend.text = element_text(size = 8))
p_b2
# ggsave("./figs/finished figs/birdphylo.tiff", p_b2, width=8, height = 8, units = "in")13.4 mammals
# subset dataset to include only mammals
mammal_data <- as.data.frame(figs_data %>%
filter(taxo_group == "mammal"))
# setting up the basic tree structure
# load tree, set node colours
mammaltree <- read.tree("./trees/mammal_species.nwk")
# prune tree to get rid of species we no longer have data for
pruned.mammaltree <- drop.tip(mammaltree, setdiff(mammaltree$tip.label, mammal_data$spp_name_phylo))
# remove underscores from tip labels
pruned.mammaltree$tip.label = gsub("_", " ", pruned.mammaltree$tip.label)
# set rownames for labelling tips
rownames(mammal_data) <- pruned.mammaltree$tip.label
# remove underscores from species name from mammal dataset
mammal_data$spp_name_phylo = gsub("_", " ", mammal_data$spp_name_phylo)
# set row names
row.names(mammal_data) <- mammal_data$spp_name_phylo
# set window
# dev.new(width=8,height=7,noRStudioGD = TRUE)
# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
# subset dataset
pers_mammal <- as.data.frame(pers_new %>%
filter(taxo_group == "mammal"))
# make this a matrix-style dataframe
pers_mammal <- data.frame(pers_mammal %>%
group_by(spp_name_phylo, personality_trait) %>%
summarise(n = n()))## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument. # remove underscores species names
pers_mammal$spp_name_phylo = gsub("_", " ", pers_mammal$spp_name_phylo)
pers_mammal <- data.frame(pers_mammal %>%
spread(personality_trait, n, fill = 0))
row.names(pers_mammal) <- pers_mammal$spp_name_phylo
pers_mammal <- pers_mammal[,2:6]
# matrix
mammal_matrix <- data.matrix(pers_mammal)
# making the tree
p_m1 <- ggtree(pruned.mammaltree, size = 0.3, layout = 'circular', branch.length = 'none') %<+% mammal_data +
geom_tippoint(aes(color=SSD_index)) +
scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") +
geom_tiplab2(size = 2.5, offset = 2, colour = "black", fontface = "italic") +
theme(legend.position = 'right')
# adding heatmap of traits
p_m2 <- gheatmap(p_m1, mammal_matrix, offset=32, width=1.3, low = "white", high = "mediumseagreen", color=NULL,
colnames=T, colnames_angle = 60, colnames_offset_y = .1, colnames_offset_x = .2) +
theme(plot.tag = element_text(size = 9, face = "bold"),
legend.text = element_text(size = 8))
p_m2
# ggsave("./figs/finished figs/mammalphylo.tiff", p_m2, width=8, height = 7, units = "in")13.5 reptiles
# subset dataset to include only reptiles
rep_data <- as.data.frame(figs_data %>%
filter(taxo_group == "reptilia"))
row.names(rep_data) <- rep_data$spp_name_phylo
# setting up the basic tree structure
# load tree, set node colours
reptree <- read.tree("./trees/reptile_species.nwk")
# prune tree to get rid of species we no longer have data for
pruned.reptree <- drop.tip(reptree, setdiff(reptree$tip.label, rep_data$spp_name_phylo))
# remove underscores from tip labels
pruned.reptree$tip.label = gsub("_", " ", pruned.reptree$tip.label)
# set rownames for labelling tips
rownames(rep_data) <- pruned.reptree$tip.label
# remove underscores from species name from mammal dataset
rep_data$spp_name_phylo = gsub("_", " ", rep_data$spp_name_phylo)
# set window size
#dev.new(width=7,height=5,noRStudioGD = TRUE)
# tree structure
p3 <- ggtree(pruned.reptree, branch.length='none', size = 0.3, layout='circular') %<+% rep_data +
geom_tippoint(aes(color=SSD_index)) +
scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") +
geom_tiplab2(align=T, linetype=NA, size=2.5, offset=4, hjust=0, colour = "black", fontface = "italic")
# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
# subset dataset
pers_rep <- as.data.frame(pers_new %>%
filter(taxo_group == "reptilia"))
# make this a matrix-style dataframe
pers_rep <- data.frame(pers_rep %>%
group_by(spp_name_phylo, personality_trait) %>%
summarise(n = n()))## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument. # remove underscores from species name from mammal dataset
pers_rep$spp_name_phylo = gsub("_", " ", pers_rep$spp_name_phylo)
pers_rep <- data.frame(pers_rep %>%
spread(personality_trait, n, fill = 0))
row.names(pers_rep) <- pers_rep$spp_name_phylo
pers_rep <- pers_rep[,2:6]
# matrix
rep_matrix <- data.matrix(pers_rep)
# add the heatmap data to our plot
rep_plot <- gheatmap(p3, rep_matrix, offset = 40, width = 3.5,
low = "white", high = "mediumseagreen", color=NULL,
colnames_position="top",
colnames_angle=60, colnames_offset_y = 0,
hjust=0, font.size=3) #just not aligning properly
rep_plot
# ggsave("./figs/finished figs/repphylo.tiff", rep_plot, width=7, height = 5, units = "in")13.6 fish
# subset dataset to include only fish
fish_data <- as.data.frame(figs_data %>%
filter(taxo_group == "fish"))
# window size
#dev.new(width=8,height=6,noRStudioGD = TRUE)
# setting up the basic tree structure
# load tree
fishtree <- read.tree("./trees/fish_species.nwk")
# prune tree to get rid of species we no longer have data for
pruned.fishtree <- drop.tip(fishtree, setdiff(fishtree$tip.label, fish_data$spp_name_phylo))
# remove underscores from tip labels
pruned.fishtree$tip.label = gsub("_", " ", pruned.fishtree$tip.label)
# set rownames for labelling tips
rownames(fish_data) <- pruned.fishtree$tip.label
# remove underscores from species name from fish dataset
fish_data$spp_name_phylo = gsub("_", " ", fish_data$spp_name_phylo)
row.names(fish_data) <- fish_data$spp_name_phylo
# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
# subset dataset
pers_fish <- as.data.frame(pers_new %>%
filter(taxo_group == "fish"))
# make this a matrix-style dataframe
pers_fish <- data.frame(pers_fish %>%
group_by(spp_name_phylo, personality_trait) %>%
summarise(n = n()))## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument. # remove underscores from tip labels
pers_fish$spp_name_phylo = gsub("_", " ", pers_fish$spp_name_phylo)
pers_fish <- data.frame(pers_fish %>%
spread(personality_trait, n, fill = 0))
row.names(pers_fish) <- pers_fish$spp_name_phylo
pers_fish <- pers_fish[,2:6]
# matrix
fish_matrix <- data.matrix(pers_fish)
# FINAL TREE
p_f1 <- ggtree(pruned.fishtree, size = 0.3, layout = 'circular', branch.length = 'none') %<+% fish_data +
xlim(-30, NA) +
geom_tippoint(aes(color=SSD_index)) +
scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") +
geom_tiplab2(size = 2.5, offset = 6, colour = "black", fontface = "italic") +
theme(legend.position = 'right')
# add the heatmap data to our plot
fish_plot2 <- gheatmap(p_f1, fish_matrix, offset = 170, width = 5.5,
low = "white", high = "mediumseagreen", color=NULL,
colnames_position="bottom",
colnames_angle=60,
hjust=0, font.size=3)
fish_plot2 
# ggsave("./figs/finished figs/fishphylo.tiff", fish_plot2, width=8, height = 6, units = "in")13.7 inverts
# subset dataset to include only inverts
invert_data <- as.data.frame(figs_data %>%
filter(taxo_group == "invertebrate"))
# setting up the basic tree structure
# load tree, set node colours
inverttree <- read.tree("./trees/invert_species.nwk")
# prune tree to get rid of species we no longer have data for
pruned.inverttree <- drop.tip(inverttree, setdiff(inverttree$tip.label, invert_data$spp_name_phylo))
# remove underscores from tip labels
pruned.inverttree$tip.label = gsub("_", " ", pruned.inverttree$tip.label)
# remove underscores from dataset and fix row names
invert_data$spp_name_phylo = gsub("_", " ", invert_data$spp_name_phylo)
row.names(invert_data) <- invert_data$spp_name_phylo
# set rownames for labelling tips
rownames(invert_data) <- pruned.inverttree$tip.label
# dev.new(width=8,height=6,noRStudioGD = TRUE)
# tree structure (cladogram, circular)
p5 <- ggtree(pruned.inverttree, branch.length='none', layout='circular') %<+% invert_data +
geom_tippoint(aes(color=SSD_index)) +
scale_color_gradient2(midpoint = 0, low = "red3", mid = "seashell2", high = "deepskyblue2") +
geom_tiplab2(align=T, linetype=NA, size=2.2, offset=2, fontface = "italic") +
theme(legend.position = "right")
# make a matrix of effect sizes (n) for each species for each personality trait to add to our plot!
# subset dataset
pers_invert <- as.data.frame(pers_new %>%
filter(taxo_group == "invertebrate"))
# make this a matrix-style dataframe
pers_invert <- data.frame(pers_invert %>%
group_by(spp_name_phylo, personality_trait) %>%
summarise(n = n()))## `summarise()` has grouped output by 'spp_name_phylo'. You can override using the `.groups` argument. # remove underscores from dataset and fix row names
pers_invert$spp_name_phylo = gsub("_", " ", pers_invert$spp_name_phylo)
pers_invert <- data.frame(pers_invert %>%
spread(personality_trait, n, fill = 0))
row.names(pers_invert) <- pers_invert$spp_name_phylo
pers_invert <- pers_invert[,2:6]
# matrix
invert_matrix <- data.matrix(pers_invert)
# add the heatmap data to our plot
invertplot <- gheatmap(p5, invert_matrix, offset = 40, width = 1.5,
low = "white", high = "mediumseagreen", color=NULL,
colnames_position="bottom",
colnames_angle=45, colnames_offset_y = 0,
hjust=0, font.size=2.5)
invertplot 
# save plot
# ggsave("./figs/finished figs/invertphylo.tiff", invertplot, width=8, height = 6, units = "in")These plots were edited together outside of R with the addition of creative commons animal silhouettes from PhyloPic to create Figures 2-6. Figure 1, the PRISMA diagram, was created using sankeymatic.com